Código
## primeras 6 líneas
head(d)Challenge #2
En esta práctica ajustaremos y validaremos un modelo de RLM, y finalmente realizaremos predicción.
The Tire Company (TTC) fabrica compuestos para llantas de F1. En la actualidad, TTC desarrolla un nuevo compuesto para la temporada 2023, con miras a superar nuevamente a su más desafiante competidor.
En el proceso de fabricación existen 3 compuestos \(x_1, x_2\) y \(x_3\) que conforman cada llanta. De acuerdo con las regulaciones, una llanta de F1 debe pesar no más de 13 kilogramos sin incluir el rin y los sensores. Además, se sabe que el costo por cada gramo de material \(j\), en dólares, es \(c_j = \sqrt{j}\), \(j=1,2,3\). TTC está convencida de que modificando la composición de la llanta es posible alcanzar un mejor rendimiento (medido en kilómetros recorridos), de tal forma que no es necesario realizar tantas paradas en pits durante una competencia regular. En la fase de experimentación, el equipo de analítica, donde usted trabaja, encontró que la formulación actual permite recorrer \(26 \pm 7\) kilómetros con un set de llantas. Se sabe que, en libras, \(x_1\in(5, 12)\), \(x_2\in(3, 20)\) y \(x_3\in(1, 6)\). Por lo tanto, el resto del peso de la llanta debe completarse con caucho sin procesar.
El día de ayer, usted fue asignado a este proyecto con miras a optimizar la composición de las llantas en TTC, de tal manera que el número de kilómetros recorridos por cada set de llantas sea considerablemente mayor que en la actualidad.
Para leer los datos hacemos:
Las primeras 6 filas de la base de datos d son
## primeras 6 líneas
head(d)Analicemos inicialmente la correlación entre las variables disponibles:
## matriz de correlación
par(mfrow = c(1,1), mar = c(.1, .1, .1, .1))
qgraph(cor(d), graph = "cor", layout = "spring",
sampleSize = nrow(d),
legend.cex = 1, alpha = 0.05)Numéricamente, la matriz de correlación es
## matriz de correlación
cor(d) y x1 x2 x3
y 1.0000000 0.29127317 -0.762485290 0.288541500
x1 0.2912732 1.00000000 0.040792570 0.049891558
x2 -0.7624853 0.04079257 1.000000000 -0.003470425
x3 0.2885415 0.04989156 -0.003470425 1.000000000
Estos resultados indican que la correlación entre \(y\) y \(x_1\) es 0.2912, entre \(y\) y \(x_2\) es -0.762, y entre \(y\) y \(x_3\) es 0.288. En cuanto a la correlaciones entre las variables independientes, podemos concluir que estas son pequeñas, lo cual sugiere que, efectivamente, \(x_1, x_2\) y \(x_3\) son independientes.
También podemos representar las correlaciones y la distribución de cada variable en un gráfico de dispersión/correlación:
## gráfico de dispersión/correlación
ggpairs(d) + theme_minimal()El gráfico 3D entre \(x_1\), \(x_2\) y \(y\) sería:
fig <- plot_ly(d,
x = ~x1,
y = ~x2,
z = ~y,
text = ~rownames(d),
color = '#BF382A')
fig <- fig %>% add_markers()
fig <- fig %>% layout(title = '\n y vs. (x1, x2)',
scene = list(xaxis = list(title = 'x1'),
yaxis = list(title = 'x2'),
zaxis = list(title = 'y')))
figEl modelo ajustado es:
## full MLR model
fit <- lm(y ~ ., data = d)
summary(fit)
Call:
lm(formula = y ~ ., data = d)
Residuals:
Min 1Q Median 3Q Max
-7.7744 -2.2090 -0.0784 1.8554 7.5166
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.24199 1.48841 16.959 < 2e-16 ***
x1 1.19280 0.15717 7.589 3.49e-12 ***
x2 -1.65137 0.08684 -19.016 < 2e-16 ***
x3 1.73492 0.26130 6.640 5.78e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.815 on 146 degrees of freedom
Multiple R-squared: 0.7584, Adjusted R-squared: 0.7535
F-statistic: 152.8 on 3 and 146 DF, p-value: < 2.2e-16
La ecuación del modelo ajustado es
\[ \hat{y} = 25.24 + 1.193 x_1 - 1.651x_2 + 1.735x_3 \]
Los intervalos de confianza del 95% para los coeficientes pueden obtenerse a través de la función confint.default() haciendo
## 95% CI para los coeficientes
confint.default(fit) 2.5 % 97.5 %
(Intercept) 22.3247667 28.159218
x1 0.8847651 1.500843
x2 -1.8215743 -1.481163
x3 1.2227821 2.247048
En términos generales, el concepto de multicolinealidad es sinómino de redundancia en las variables independientes.
Uno de los supuestos fuertes del modelo de RLM es que las variables \(X_1, X_2,\ldots,X_n\) son independientes. Cuando esto no ocurre, los estimadores de \({\beta} = (\beta_1,\beta_2,\ldots,\beta_k)\) tienen propiedades distintas a las ya conocidas.
Desde el punto de vista formal, la existencia de multicolinealidad puede probarse utilizando el ill-condition number (ICN)
## ICN
kappa(fit) [1] 44.46287
Teniendo en cuenta que el ICN es \(>30\), aparentemente, existe multicolinealidad entre \(x_1, x_2\) y \(x_3\).
Ahora, si estamos interesados en determinar cuál de la(s) variable(s) independiente(s) con mayor grado de colinealidad, utilizamos el VIF:
## VIF
car::vif(fit) x1 x2 x3
1.004185 1.001697 1.002526
Al analizar el VIF, es claro que ningún valor es \(>5\).
Recientemente, se han implementado otras pruebas de multicolinealidad en el paquete mctest:
## otras pruebas de multicolinealidad
mctest(fit)$odiags results detection
Determinant 0.99582063 0
Farrar Chi-Square 0.61635301 0
Red Indicator 0.03726144 0
sum of Lambda Invers 3.00840801 0
Theil Indicator -1.50847225 0
Condition Number 15.61480762 0
De acuerdo con estos resultados, no es posible hablar de la existencia de multicolinealidad. En parrticular, no aparece el valor de 1 en la columna detection.
La validación de los supuestos puede hacerse via valores \(p\) como se muestra a continuación.
Normalidad
## prueba de Normalidad
shapiro.test(rstudent(fit))$p.value[1] 0.8899012
Como el valor \(p\) es \(>0.05\), entonces los errores del modelo siguen una distribución Normal.
Varianza constante
## prueba de varianza constante
car:::ncvTest(fit)$p[1] 0.3356801
Como el valor \(p\) es \(>0.05\), podemos concluir que los errores tienen varianza constante.
Independencia
## prueba de independencia
car:::durbinWatsonTest(fit)$p[1] 0.302
Este resultado indica que los errores del modelo ajustado son independientes.
Para identificar outliers usamos los residuales estudentizados del modelo:
## outliers?
r <- rstudent(fit)
res <- which(r > 3 | r < -3)
nout <- ifelse(length(res) == 0, 0, length(res))Este resultado indica que hay 0 outliers en los datos.
Para identificar este tipo de observaciones, utilizamos la distancia de Cook.
En R procedemos de la siguiente manera:
## gráfico de la distancia de Cook
plot(fit, which = 4, las = 1)Este resultado indica que las observaciones 44, 57 y 142 podrían considerarse influenciales.
Ahora, podemos encapsular las pruebas anteriores en una función:
## cargar función
source('https://www.dropbox.com/scl/fi/1z0gsv6g4a9eowcoqs4xf/validar_supuestos.R?rlkey=ioqrt3r3sl3vh7wa7sc82gjqh&dl=1')Esta función fue creada por el profesor Jorge Vélez, PhD
<https://jorgeivanvelez.netlify.app/> para la validación
de los supuestos básicos en RLM.
En caso de presentar algún inconveniente, por favor
repórtelo a <jvelezv@uninorte.edu.co>
Última modificación: Abril 15, 2024
Finalmente, utilizamos la función validar_supuestos() para validar todos los supuestos a la vez e identificar outliers:
## validar supuestos
## fit es el objeto que contiene el modelo de RLM
validar_supuestos(fit) outliers normalidad homocedasticidad independencia
"no hay outliers" "cumple" "cumple" "cumple"
Conclusión: Los supuestos de Normalidad, independencia y varianza constante de los errores se cumplen. Por lo tanto, el modelo es válido para predecir. Además, parecen no existir outliers en los datos.
validar_supuestos()A continuación, cargaremos la función que hemos creado en nuestro repositorio de GitHub, la cual es una herramienta más completa para la validación y construcción de un modelo de regresión lineal múltiple. En ella seguiremos una serie de pasos para la construcción del modelo, omitiendo el primer paso de ajuste inicial.
Para utilizar esta función en nuestro entorno de trabajo, primero cargaremos el código desde el repositorio:
source("https://raw.githubusercontent.com/unfresh25/Challenge--2/master/function.R")Esta función fue creada por Jorge Borja
<https://jorgeborja-portfolio.vercel.app/> para la validación
de la construcción de un modelo de regresión lineal múltiple.
En caso de presentar algún inconveniente, por favor
repórtelo a <jborja@uninorte.edu.co>
Última modificación: Junio 05, 2024
Luego, procederemos a evaluar el modelo (fit) utilizando la función creada, en este caso utilizamos plot = T para ver el gráfico de la distancia de cooks. Entonces,
validar_supuestos(fit, plot = T)1. Validación
Significancia Global Significancia Marginal R2 ajustado
"Significativo" "Significativas" "Significativo"
2. Verificación de residuales
Normalidad Independencia Homocedasticidad
"Cumple" "Cumple" "Cumple"
3. Multicolinealidad
Multicolinealidad
"Cumple"
4. Outliers/Influenciales
Outliers Influenciales
"No hay outliers" "44, 57, 142"
Si fuese de interés determinar el rendimiento de una llanta con compuestos
\[ \mathbf{x}_0= (10, 4, 6) \]
procedemos de la siguiente forma:
## x0
x0 <- data.frame(x1 = 10, x2 = 4, x3 = 6)
x0## estimación
predict(fit, newdat = x0) 1
40.97405
Por lo tanto, \(\widehat{E[Y|\mathbf{x}_0]} = 40.97\).
interval a predict():El intervalo de confianza del 95% es
## confidence interval
predict(fit, newdat = x0, interval = 'confidence') fit lwr upr
1 40.97405 38.88432 43.06377
Así las cosas, el promedio poblacional de la distancia recorrida cuando se utiliza la mezcla de componentes \(\mathbf{x}_0= (10, 4, 6)\) es
\[ \mu | \mathbf{x}_0 \in (38.88, 43.06) \] con una confianza del 95%.
Finalmente, el intervalo de predicción del 95% es
## prediction interval
predict(fit, newdat = x0, interval = 'prediction') fit lwr upr
1 40.97405 35.03113 46.91697
Por lo tanto, la distancia recorrida por el próximo set de llantas en el que se utilice la mezcla de componentes \(\mathbf{x}_0= (10, 4, 6)\) será
\[ Y | \mathbf{x}_0 \in (35.03, 46.92) \]
Para más información sobre el cálculo de los intervalos de confianza y predicción a partir de un modelo de RLM ajustado, se recomienda consultar ?predict.lm.
Para más detalles sobre el modelo de RLM se sugiere consultar el Capítulo 3 del texto Modelos de Regresión: Una aproximación práctica con R.